% The code is based on the following paper:
%   Chen, T. Y., Lin, Y. L., and Tzeng, L. Y., forthcoming.
%   Estimating probability weighting functions through option pricing bounds. Review of Asset Pricing Studies.
%
% Copyright: Tzu-Ying Chen, Yo-Lan Lin, Larry Y. Tzeng
% Date: January 30, 2024

clear

%% Setting
AllYear = repmat(1981:2022, 12, 1);                        
AllMonth = repmat([1:12]', 1, 2022 - 1981 + 1);            
AllDate_EXP = nweekdate(3, 6, AllYear(:), AllMonth(:));                                             
Target_AllDate = [str2num(datestr(AllDate_EXP - 30, 'yyyymmdd')) ...
                  str2num(datestr(AllDate_EXP, 'yyyymmdd'))];

Index_Drop = find((fix(Target_AllDate(:, 1) / 10000) < 1981) | ...
                  (fix(Target_AllDate(:, 1) / 10000) > 2021));
Target_AllDate(Index_Drop, :) = [];                                        
clear Index_Drop AllYear AllMonth AllDate_EXP Target_TTM

%% Load Data
FileName = ['StockReturn_19262021_SP500.txt'];
Data_RET = load(FileName);    
Data_RET = sortrows(Data_RET, 1);                                          
clear FileName

RET_SP500 = nan(size(Target_AllDate, 1), 1);                               
for d = 1:size(Target_AllDate, 1)
    Target_Date = Target_AllDate(d, 1);    
    Index_Begin = find(Data_RET(:, 1) <= Target_Date);       
    Index_Begin = Index_Begin(end);                                        

    Target_ExDate = Target_AllDate(d, 2); 
    Index_End = find(Data_RET(:, 1) < Target_ExDate);                         
    Index_End = Index_End(end);                                             

    RET_SP500(d) = prod(1 + Data_RET((Index_Begin + 1):Index_End, end)) - 1;
    clear Index_Begin Index_End Target_Date Target_ExDate
end
clear Target_AllDate Data_RET
clear d

%% Calculate Probability from Continuous to Discrete State Space (Constantinides, Jackwerth, and Perrakisy, 2009)
Data = 1 + RET_SP500;                                               
h = power(4 / 3, 1 / 5) * std(Data) * power(length(Data), - 1 / 5);        
  
State = (- 0.6):0.02:0.6;                                                  
State = exp(State);                                                        
State_PROB = nan(size(State));                                            
for n = 1:length(State)
    if (n > 1) & (n < length(State))
        Point_Begin = interp1([0 1], ...
                              [State(n - 1) State(n)], ...
                              0.5);
        Point_End = interp1([0 1], ...
                            [State(n) State(n + 1)], ...
                            0.5);
    elseif n==1
        Point_End = interp1([0 1], ...
                            [State(n) State(n + 1)], ...
                            0.5);
        Point_Begin = State(n) - ...
                      (Point_End - State(n));
    elseif n==length(State)
        Point_Begin = interp1([0 1], ...
                              [State(n - 1) State(n)], ...
                              0.5);
        Point_End = State(n) + ...
                    (State(n) - Point_Begin);
    else
    end

    Smooth_Data = linspace(Point_Begin, Point_End, 100000);
    Smooth_Density = ksdensity(Data, Smooth_Data, ...
                               'kernel', 'normal', ...
                               'bandwidth', h, ...
                               'function', 'pdf');    
    State_PROB(n) = trapz(Smooth_Data, Smooth_Density);
    clear Point_Begin Point_End Smooth_Data Smooth_Density
end
State_PROB = State_PROB / sum(State_PROB);
clear n

Index_Drop = logical(State_PROB < 0.00001);
State(Index_Drop) = [];                                                    
State_PROB(Index_Drop) = [];                                               
State_PROB = State_PROB / sum(State_PROB);                                 
clear Index_Drop Data h           

%% De-Mean and Reintroduce Risk Premium Equal to 4% (Annualized) in Each Period
FileName = ['AllTradingDate19962021_Monthly.mat'];
load(FileName, ...
     'Target_AllDate');        
clear FileName

Date_Monthly = Target_AllDate;
State_Monthly = nan(size(Date_Monthly, 1), length(State));                 
for d = 1:size(Date_Monthly, 1)
    Target_Date = Date_Monthly(d, 1);
    Target_TTM = Date_Monthly(d, 3);

    FileName = ['OP' num2str(Target_Date) '_TTM' num2str(Target_TTM) '.mat'];
    load(FileName, ...
         'Data', 'Index_TTM', 'Index_RF')
    clear FileName    

    TTM = Data(1, Index_TTM) / 365;                                        
    RF = Data(1, Index_RF);                                                
    State_Monthly(d, :) = State - sum(State .* State_PROB) + ...
                          power(1 + 0.04 + (exp(RF) - 1), TTM);
    clear Target_Date Target_TTM Data Index_TTM Index_RF                
    clear TTM RF 
end
clear Target_AllDate
clear d

%% Output
FileName = ['GrossRET_1981_2021_DiscreteState'];
save(FileName)
clear FileName
